analyticalVar <- function(X, Y, n, n_omega, beta_tilde, S) {
          H_R_l <- lapply(1:nrow(X), FUN = function(i) {
            delta <- X[i, ]*as.numeric((Y[i] - t(X[i, ])%*%beta_tilde)) + S%*%as.matrix(beta_tilde)
            return((delta%*%t(delta)))
          })
          H_R <- apply(array(unlist(H_R_l) , c(ncol(X),ncol(X), nrow(X)) ), 1:2, sum)
          var_robust <- solve(n_omega)%*%H_R%*%solve(n_omega)
          se_robust <- sqrt(diag(var_robust))
    return(se_robust)
}

timeVar <- function(N, S = 2) {
  
  Y_err <- 2
  corr <- 2
  
  # Draw the true X's
  Z1_mean <- 7
  Z2_mean <- 9
  Z1 <- rpois(N, Z1_mean)
  Z2 <- Z1*corr + rpois(N, Z2_mean)
  b0 <- 10
  b1 <- 12
  b2 <- -3
  Y_res <- rnorm(N, mean = 0, sd = Y_err)
  Y <- b0 + b1*Z1 + b2*Z2 +  Y_res
  V1 <- rnorm(N, 0, 1)
  V2 <- rnorm(N, 0, 1)
  X1 <-  Z1 + V1
  X2 <- Z2 + V2
  X <- cbind(1, X1, X2)
  
  x_prime_x <- t(X)%*%X 
  x_y <- t(X)%*%Y
  beta_tilde <- c(b0, b1, b2)
  S <- matrix(0, 3, 3)
  diag(S) <- c(0, 1, 1)
  n_omega <- t(X)%*%X - N*S

  
  sim_time <- system.time(code <- varianceMVN(X = X, Y = Y, sigma_sq = Y_err, nsims  = 500, coef = beta_tilde, S, N, x_prime_x, x_y))

  analytical_time <- system.time(code <- analyticalVar(X, Y, N, n_omega, beta_tilde, S))

  return(data.frame(simulation = sim_time['elapsed'], analytical = analytical_time['elapsed'], N))
}

N <- c(100000, 500000, 1000000, 2000000, 3500000, 5000000)


time_variance <- as.data.frame(do.call(rbind, lapply(N, FUN = function(n) timeVar(n))))

save(time_variance, file = "simulation_output/time_test_analytical.Rdata")

